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Synopsis 

The behavior of suspensions of rigid particles in a non-Newtonian fluid is studied in the 
framework of a nonlinear homogenization method. Estimates for the overall properties of 
the composite material are obtained. In the case of a Herschel-Bulkley suspending fluid, it is 
shown that the properties of a suspension with overall isotropy can be satisfactory modeled 
as that of a Herschel-Bulkley fluid with an exponent equal to that of the suspending fluid. 
Estimates for the yield stress and the consistency at large strain rate levels are proposed. 
These estimates compare well to both experimental data obtained by Mahaut et al. (2007) 
and to experimental data found in the literature. 

I Introduction 

Heterogeneous systems consisting of particles suspended in a fluid medium make up in a wide 
variety of materials of practical interest, both natural (slurries, debris flows, lavas, ...) or man- 
made (concretes, food pastes, paints, cosmetics, ...). This abundance explains why the behavior 
of these materials have been extensively studied from both a theoretical and an experimental 
point of view. 

In the framework of man-made materials, it does not seem to exist a well established method 
to obtain a material with given rheological characteristics from components having known prop- 
erties despite the fact that is is a problem encountered in numerous industrial processes. The 
set-up of such a methodology requires the ability to predict the overall behavior of the material 
from that of its constituents. A general answer for this problem has yet to be found. 

The main difficulty in modeling the behavior of suspensions comes from the fact that the 
material is multiscale and contains many interacting constituents. Experimental data show that 
the overall rheological properties of a suspension depend upon the shape and size of the particles, 
the interaction between particles (colloidal or noncolloidal), the interaction between the particles 
and the suspending fluid (hydrodynamic) , the properties of the suspending fluid (Newtonian or 
not) and the type of flow the suspension is subjected to. 

In this paper, we focus on suspensions made up of noncolloidal particles dispersed in a yield 
stress fluid such as suspensions made up of a coarse fraction and a colloidal fraction. If the 
colloidal particles are much smaller than the coarse ones, the latter interact with the other 
components of the suspension only through hydrodynamic interactions: they see a homogeneous 
phase fluid which behavior is that of the colloidal suspension [Sengun and Probstein ( 1989a, b); 
Ancey and Jorrot (2001)]. 

Like for any other suspensions, the rheological properties of non-Newtonian suspensions de- 
pend upon the shape, the surface texture and the size distribution of the coarse particles. 
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In addition to experimental approaches, modeling the behavior of non-Newtonian suspen- 
sions using their constituents characteristics in the framework of numerical methods appears 
promising because of the ability of these methods to track the particles localization in the sus- 
pension and to account for several interparticle interactions. Solving the problem for noncolloidal 
particles immersed in a Newtonian fluid is feasible using modern computers and provides useful 
results [Brady (2001)]. The most serious drawback of these methods to their generalization to 
non-Newtonian suspending fluid is the fact that the Newtonian hydrodynamic force applied to 
the particle is evaluated by means of close form equations depending on the macroscopic strain 
rate and on the particle velocity relative to the bulk fluid [Bossis and Brady (1984); Brady and 
Bossis (1988)]. Thanks to these relations, it is not necessary to solve the continuous Stokes 
equation in the bulk fluid in order to simulate the behavior of the suspension. This enables 
to simulate the flow of representative elementary volume of a suspension containing numerous 
particles at a reasonable computational cost. To our knowledge, such an equation does not exist 
for particles immersed in a yield stress fluid, which prevents a generalization of the method to 
be easily proposed. 

Trying to evaluate the velocity field in the bulk yield stress fluid by means of a finite element 
method is such a time consuming task that no situation of practical interest can be solved using 
these tools today [Johnson and Tezduyar (1997); Roquet and Saramito (2003); Yu and Wachs 
(2007)]. 

To obtain estimates for the overall characteristics of non-Newtonian suspensions, the change 
of scale method appears then to be a powerful tool. As a reminder homogenization technique 
aims at identifying the macroscopic properties of a material modeled as a continuous medium 
from those of their constituents. The first results in this field were obtained by Einstein (1906) 
for the viscosity of a dilute suspension. Since, different problems have been studied such as 
the viscosity of multimodal suspensions [Farris (1968)], the viscosity of concentrated suspensions 
[Krieger and Dougherty (1959); Frankel and Acrivos (1967)], the effect of the shape and orien- 
tation of particles on the behavior of Newtonian suspensions [Batchelor (1971)], the effect of 
interparticle interactions [Batchelor and Green (1972)] or the effect of Brownian motions on the 
overall behavior [Batchelor (1977); Russel et al. (1995)]. 

It is worth noticing that it has not been possible to obtain exact solutions for most of the 
problems cited above; generally, only estimates of the overall characteristics of the suspension 
have been obtained. This situation is very similar to the one prevailing in the field of homog- 
enization approaches to the behavior of solid materials [Zaoui (2002)]. This is not surprising 
since the two problems are very similar as it has been recognized by Batchelor and Green (1972). 
The main difference comes from the fact that the morphology of the heterogeneities within the 
representative elementary volume is a given for solid materials whereas it is an unknown for 
particle suspensions since the fiow of the representative elementary volume of suspension and 
the morphology of the particles are coupled 

Interestingly, novel developments have taken place in nonlinear continuum micromechanics 
in the last twenty years so that both estimates and variational bounds are now available [see 
Suquet (1998) for a review]. Because of the similarities between the solid problem with the liquid 
one, it is quite natural to address the modeling of the behavior of non-Newtonian suspensions 
within this framework. 

The aim of this paper is to provide a first approach to the overall behavior of a suspension 
of non-Brownian and noncolloidal particles immersed in an incompressible yield stress fluid. 
First, the main features of the homogenization approach to the behavior of non-Newtonian 
suspensions are recalled. Then, the secant-method of Castaheda (1991) and Suquet (1993) is 
applied in order to estimate the overall behavior of the suspension. The estimates are compared 
to new experimental results obtained by Mahaut et al. (2007) and to experimental results of the 
literature in the next part. Finally, the validity of the theoretical approach is discussed before 
some general conclusions are drawn. 
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II Homogenized behavior 



We examine the overall behavior of a suspension, the liquid phase of which is homogeneous and 
nonlinear. We restrict ourselves to the situation where the constitutive behavior of the fluid is 
characterized by an energy function w (d) where d = \J'ld : d denotes the second invariant of the 
Eulerian strain rate tensor d. The Cauchy stress cr is obtained by differentiation of the potential 
w with respect to the strain rate tensor d if tt) is differentiable. Otherwise, the derivative should 
be interpreted as the subdifferential of convex analysis. 

The condition of naught strain for the particles can also be written in term of a dissipation 
potential with w defined by 

w; (d) = if d = -u; (d) = oo if d / (1) 
Thus, the suspension is made of a heterogeneous medium, the behavior of which is described 

by 

.. = ||«i,s) (2) 

w being a strict convex function. The location in the representative elementary volume is defined 
by the position vector x_. 




u(x) = D ■ X 



n 



Figure 1: The representative elementary volume of the suspension submitted to a macroscopic 
strain rate loading (Hashin boundary condition). 



It is assumed that it is possible to define a representative elementary volume of the suspension 
occupying a domain Q with boundary dQ such that it is large enough to be of typical composition 
and its overall properties do not depend on the way it is loaded at the macroscopic scale, fi^ 
and Qi denote the solid and the liquid domain respectively. The volume fraction of particles ip 
is the ratio of the volume fraction of $7^ in 0. For simplicity, it is assumed that the boundary of 
0, is located in the liquid domain as depicted in Fig. 1. At the microscopic scale, Q is considered 
as a structure. < a > (resp. < a >a with a = s, i) denotes the average of a over Q (resp. f^a). 
The liquid phase is homogeneous. We adopt an Eulerian description of the movement and we 
restrict our attention to the situations where the evolutions of the system are quasistatics (i.e. 
inertial effects are negligible) and all the long range forces other than the hydrodynamic ones 
are negligible. 

As shown by Hill (1963), the overall behavior of the suspension reads: 

dW 

Y: = ^—(D) with W (D) = min (w (d)) (3) 

where C {D) denotes the set of Eulerian strain rate fields kinematically admissible with D. It is 
recalled that a strain rate field d defined over is said to be kinematically admissible with D if 
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exists a velocity field u defined over 0, complying with the Hashin condition 

(yxGdn) u{x) = D-x (4) 

such as: 

1 ^ ' - t„ 



(V X G il) d{x) = - l^grad u (x) + grad u {x)j (5) 

When the fluid obeys a Herschel-Bulkley law with a yield stress Tc, a consistency r] and a 
power law exponent n > 0, the dissipation potential reads: 



w (d) = Tcd+ -^(t^^ if trd = 
^ ' n + 1 

w (d) = oo if trd / 

Then the fluid's state equation reads 



(6) 



d = if y/s : s/2 < r, 

s= (^Tc + r]d"-^ i trd = if a/s : s/2 > t, 



(7) 



where p = — trcr/3 is the hydrostatic pressure, s = a+p5 the deviator of a and 6 the second order 
unit tensor. As it is common for incompressible materials, the pressure p is not determined by 
the state law. Considering that the particles are rigid and that the bearing fluid is homogeneous 
and incompressible, it is easily shown from Eq. (3) that the macroscopic potential is also defined 
by 



W (D) = min (1 - w) 
deg{D) 



d) +^(d-^') 

n + 1 \ 11 



if trU = (8) 



where Q (D) denotes the subset of C (D) which elements comply with the naught strain rate 
constraint over the domain occupied by the particles and the incompressible constraint over the 
fluid domain. Of course, the set Q (D) is only deflned for macroscopic strain rate complying with 
the condition tiD = 0. Eq. (8) is completed by the condition W (-D) = oo if tiD ^ 0, which 
enforced the incompressible constraint at the macroscopic level. 

The set Q (D) being convex, the minimization problem (8) admits only one solution, which 
ensures the validity of the method. The identification of the macroscopic behavior of the suspen- 
sion from Eq. (3) or Eq. (8) requires the resolution of a continuous convex minimization problem. 
This problem of minimization has to be solved for each morphology of the suspension defined by 
the shape of the particles and the distribution of the particles within the suspension. Of course, 
it is not possible to solve this problem in most situations of practical interest. 

To remedy this difficulty, various estimation techniques of the macroscopic behavior have 
been proposed, in particular, by Castaheda (1991, 1996, 2003) and Suquet (1993). 

The key feature in these methods is the use of a rigorous variational principle (Eq. 3 for 
the yield stress fiuid suspensions) to determine the best possible choice of a "linear comparison 
composite" to estimate the effective behavior of the nonlinear one. A detailed description of these 
methods is beyond the scope of the present paper. The reader is referred to [Suquet (1998)] for 
a more detailed review. 

In the following, we obtained estimates relevant to our problem in the framework of a simple 
approach which does make use explicitly of the variational Eq. 3. The main features of the 
method used to obtain these estimates are recalled in the following section (largely inspired by 
the presentation of Suquet (1997)) for completeness of the paper. 
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Ill Secant estimate of the behavior 



It is assumed that the sohd particles are isotropicaUy distributed over the representative elemen- 
tary volume. The macroscopic behavior of the suspension is therefore also isotropic. 
It is possible to write the behavior of the nonlinear fluid in the following form 



(9) 



where denotes the secant modulus of the fluid phase. The secant modulus is no more 

than the apparent viscosity of the fluid, as depicted in Fig. 2. 



Sij{d) 




Figure 2: Secant modulus formulation of nonlinear incompressible isotropic fluid. 



Using Eq. (9) for the state equation of the suspending fluid, it is possible to replace the 
original nonlinear problem of homogenization by a linear homogenization problem for a suspen- 
sion of rigid particles immersed in a heterogeneous incompressible isotropic fluid with apparent 



viscosity r] (x) = fi^^^ yx, d [x 

Then, the difficulty which remains to be solved is to calculate, or to estimate, the field d 
over the domain occupied by the fluid phase in the representative elementary volume. Since it 
is impossible to analytically determine the local response of the nonlinear fluid, it is impossible 
to compute the secant modulus field over the representative elementary volume. Then an ap- 
proximation has to be introduced to make analytic calculations feasible. The approximation for 
secant methods consists of replacing the secant modulus field over the representative elementary 
volume by a modulus which is uniform over subdomains of the representative elementary volume. 
In this paper, we consider only one domain for the liquid phase to simplify. Then the secant 
modulus is uniform over the fluid phase. This estimate reads: 



(V X G f]^ ) (x,d{x)]^ (df^^ 



(10) 



where df^ is an equivalent effective strain rate which remains to be defined as a function of the 
mean value of the field d in the fluid phase, and thus, of the value of the macroscopic strain 
rate. The replacement of the field of heterogeneous secant modulus defined over the fluid phase 
by a homogeneous fleld simplifies considerably the resolution of the problem. Consequently, the 
estimate technique of the nonlinear macroscopic behavior includes three steps: 

1. First of all, it is necessary to solve a linear homogenization problem for a suspension of 
particles immersed in an isotropic homogeneous fiuid of viscosity fie. If the suspension is 
isotropic at the macroscopic scale, the macroscopic behavior is characterized by a macro- 
scopic viscosity proportional to the viscosity of the fiuid, the coefficient of proportionality 
depending on the morphology of the particles. As this problem has been the subject of 
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many works, numerous results and estimates are available in the literature dealing with 
the rheology of suspensions or the homogenization to the behavior of heterogeneous lin- 
ear materials. Let g{(p) denotes the coefficient of proportionality linking the macroscopic 
viscosity to the microscopic viscosity. 

fi^'^^ = fieg{^) (11) 

2. Then, one must choose a measure of the effective strain rate for the fluid phase. As the 
considered material is isotropic at the microscopic scale, the strain rate is characterized by 
the second order moment of the quantity d defined by: 



df=\/{d'h (12) 

This choice corresponds to the modified approach described by Suquet (1997). Castafieda 
(1991) demonstrated that this choice is optimal in the framework of a variational approach 
to the solution of the nonlinear homogenization problem under consideration. It would 
have been simpler to choose df^ = < d >£, a quantity easily computed from the egality 
(d)^ = D/{l — {p). Unfortunately, Suquet (1997) showed that estimates of the overall 
properties of the heterogeneous material obtained using this effective liquid strain rate are 
less accurate than those obtained using Eq. (12). 

3. Finally, the nonlinear character of the problem is taken into account by integrating into 
the relation of homogenization (11) linking to the viscosity of the liquid phase, the 

fact that the value of the fluid viscosity depends on df^. As df^ depends on the value of 
the macroscopic strain rate D, the value of the macroscopic secant modulus also depends 
on the value of D. The only difficulty implementing this step is the calculation of df^ as a 
function of D for the particular homogenization scheme used. It has been shown by Kreher 
(1990) that: 

{d')e = -^^D' = -^g{^)& (13) 

1 — if OfJ.i 1 — if 

By remplacing fi£ by fi^^^ in the relation of linear homogenization Eq. (11) and then by combining 
the obtained equation with the localization Eq. (13), one obtains the following estimate for 
the macroscopic secant modulus (i.e. apparent viscosity) of any nonlinear materials with any 
isotropic microstructure: 



H^""" (if, d) = g{^) * fi"""^ (d) with d = dJ^^ (14) 



1-ip 

Using notations of Sengun and Probstein (1989b), the apparent viscosity Eqs. (14) reads 



r]{ip, 7) = -qcri^p) * Vlriies) with 7eff = ji ^^^^ (15) 

V ^-f 

which is much more general than Eqs (2.8), (2.9), (4.1) and (4.2) of Sengun and Probstein (1989b) 
as they allow to take into account any estimate g{f) {i.e. r]cr{^)) of the relative viscosity of a 
Newtonian suspension. It is worth noting that Eqs. (14) and (15) are valid for any particles 
shape and dispersity. Interestingly, even if the estimate ^(v?) does not rely on a morphological 
model, Eq. (13) allows to estimate the localization factor associated with the relative viscosity 
function ^((/j). 

Coming back to the Herschel-Bulkley suspension problem, the secant modulus of the sus- 
pending fluid reads: 

[x,d{x)) = ^ + {d{x)Y^' (16) 
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Then, putting Eq. (16) into Eq. (14) yields 



D 

with: 



—horn 

/i^'""^ ((/?, D) = + n^^^D""-^ (17) 



rr- = r,V(l - ^)g{^) (18) 



^hom 
' c 



and: 

T.-1 

2 



(19) 



[1 - 99_ 

It is therefore predicted that, at the macroscopic scale, the suspension behaves as a Herschel- 
Bulkley fluid with same exponent as that of the suspending fluid. Overall yield stress and 
macroscopic consistency are defined by Eqs. (18) and (19). This result does not depend on the 
scheme (i.e. g{(p)) used to link the viscosity of the bearing fluid to the overall viscosity of the 
suspension. The quality of the prediction depends thus only on the validity of the assumption 
set I ^ J /^\\ „„„ u„ u,, ^, ,,sct/3eS\ 



that the field yx,d{x)j can be estimated by the quantity fi^ {d 

IV Experimental validation 

An experimental procedure has been designed by Mahaut et al. (2007) which complies with 
the assumptions made to obtain the theoretical results presented above. As this procedure 
is described in detail in the paper entitled Yield stress and elastic modulus of suspensions of 
noncolloidal particles in yield stress fluids, no further details concerning the experimental work 
are given in this paper; we restrict ourselves to comparisons between experimental data and 
theoretical predictions. 



A Elastic modulus vs. yield stress 

The accuracy of the estimates obtained in the framework of the theoretical approach presented 
above for the overall properties of the yield stress suspension depends on the assumption made 
to take into account the nonlinear behavior of the suspending fiuid and on the scheme used to 
estimate the overall linear behavior of the suspension. It is possible to check experimentally the 
validity of the assumption made on the heterogeneities of the secant modulus over the liquid 
domain irrespective of the errors induced by the choice of a particular homogenization scheme. 
For this, it is enough to remark that Eq. (11) enables to calculate the macroscopic elastic modulus 
^hom q£ g^j^ isotropic suspension of particles dispersed in an isotropic incompressible linear elastic 
matrix whose shear modulus is equal to G (both problems pose exactly in the same way provided 
that d and /x^ be identified with the infinitesimal strain tensor and the elastic shear modulus). 
As a consequence, it is possible to obtain a general relationship between the dimensionless elastic 
modulus and the dimensionless yield stress of a suspension of rigid particles dispersed in a yield 
stress fluid that is true whatever the scheme as long as the particle distribution is isotropic and 
using an uniform secant modulus estimate is relevant. Combining Eqs. (11), (18) and (19) yields 
the relations: 

r^^/r, = ^(l-(^)Gh°-/G (20) 

and 

Moreover, the yield stress and the consistency not being independent of one another, it is also 
possible to determine the consistency from the yield stress and the concentration: 

^hon./^ = ^^'^ '^'^ (22) 
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In Fig. 3, we have plotted the dimensionless yield stress TcItc as a function of the dimensionless 
quantity \/{\- ^p)G'^°'^/G for all the systems studied by Mahaut et al. (2007). It is recalled 
that yield stress fluids have a solid linear viscoelastic behavior below the yield stress, so that 
the macroscopic elastic modulus of the suspensions could be experimentally measured through 
oscillatory shear measurements. A good agreement between the experimental results and the 




_2 1 1 1 1 1 

-2 -1.5 -1 -0.5 0.5 1 

log - ^)G'hom/G! - l) 

Figure 3: Dimensionless yield stress r^°™/rc as a function of {1 — ip)G^°^^ / G for all the systems 
studied by Mahaut et al. (2007). The open square symbols are suspensions of polystyrene and 
glass beads in bentonite, the open circle symbols are suspensions of glass beads in carbopol while 
the solid diamond symbols are suspensions of polystyrene and glass beads in emulsion. The 
figure's coordinates were chosen so that the y = x line represents the theoretical relation (20). 

micromechanical estimation (20) (which is plotted as a straight line y = x in these coordinates) 
is observed. These results show that the data are consistent with the assumption that an uniform 
estimate of the secant modulus over the fluid domain allows to accurately estimate the overall 
properties of the suspension in the studied situations. 



B Elastic modulus 

In this section, we summarize the results of the elastic modulus measurements performed on all 
the materials. 

The evolution of the dimensionless modulus as a function of the volume fraction if of 

noncolloidal particles for all the studied materials are summarized in Fig. 4. It is observed that 
the experimental data are very well fitted to the Krieger- Dougherty law (Krieger and Dougherty 
(1959)): 

^hom / ^ — 2.5i/?r, 



g{^) = 1 - (23) 

G \ ipmj 

The value of the maximum packing fraction = 0.57 was fixed by means of a least squares 
method. This value is very close to the value ipm = 0.605 measured locally very recently in 
dense suspensions of noncolloidal particles in Newtonian fiuids by Ovarlez et al. (2006) through 
MRI techniques. The small discrepancy between the two values of tpm comes certainly from the 
anisotropy induced by the fiow in the experiments of Ovarlez et al. (2006) as it was also observed 
in the experiments performed by Parsi and Gadala-Maria (1987). 



C Yield stress 

We now present the results of the yield stress measurements. 
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Figure 4: Dimensionless elastic modulus vs the beads volume fraction 99 for all the 

systems studied by Mahaut et al. (2007). The open square symbols are suspensions of polystyrene 
and glass beads in bentonite, the open circle symbols are suspensions of glass beads in carbopol 
while the solid diamond symbols are suspensions of polystyrene and glass beads in emulsion. 
The solid line is the Krieger-Dougherty Eq. (23) with = 0.57. 



The yield stress Tc was measured with a method which avoids destroying the homogeneity 
and the isotropy of the suspension. The evolution of the dimensionless yield stress T^°™/rc as 
a function of the solid volume fraction of noncolloidal particles is depicted in Fig. 5 for the 
studied materials. 

Fig. 5 clearly shows that the yield stress of the suspensions reads as the product of the 
suspending fluid yield stress times a function /(y?). When both the dimensionless yield stresses 
and the dimensionless elastic moduli are drawn on the same diagram (see Fig. 6), it is obvious 
that the yield stress function /(</?) is different from the elastic modulus function g{'~p). 
It is possible to directly use Eqs. (17), (18) and (19) to evaluate the function /((/j)- 
For the dilute suspensions, the Einstein relation g^^{(p) = l + 5/2(/? is exact to the first order 
of ip. Introducing this relation in Eqs. (18) and (19) and keeping only the first order terms yield 
the dilute estimates: 

= re(l + 3/4(^) (24) 



The coefficients of growth of these two estimates with the solid volume fraction are different 
from that of the Einstein law. According to the experimental results depicted in Fig. 6 the 
change of scale method also predicts that the two functions g{ip) and f{ip) are different. The 
only approximation performed to obtain Eqs. (24) and (25) is to assume that a uniform estimate 
of the secant modulus over the fluid domain enables to accurately estimate the overall properties 
of the suspension. Thus, estimates (24) and (25) are not rigorously exact in contrast with the 
Einstein law. It is worth noticing that for a suspension of particles in a Bingham fluid, the overall 
consistency of the suspension is given by the classical Einstein function 1 + 5/2^9. 

For larger values of the solid volume fraction, the viscosity of a Newtonian suspension is classi- 
cally estimated using the Krieger-Dougherty equation [Krieger and Dougherty (1959); Quemada 
(1985)]. Putting the second equality of Eq. (23) into Eqs. (18) and (19) yields the Krieger- 
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Figure 5: Dimensionless yield stress r^°™/rc vs the beads volume fraction ip for all the systems 
studied by Mahaut et al. (2007). The open square symbols are suspensions of polystyrene and 
glass beads in bentonite, the open circle symbols are suspensions of glass beads in carbopol 
while the solid diamond symbols are suspensions of polystyrene and glass beads in emulsion. 
The solid line is the theoretical prediction (26) with </?^ = 0.57. The dashed curve is the dilute 
estimate (24) for the yield stress. 
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Figure 6: Dimensionless yield stress r^°™/rc (empty symbols) and dimensionless elastic modulus 
QhoTci jQ (^golid symbols) vs the beads volume fraction Lp for all the suspensions studied by Mahaut 
et al. (2007). 
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Dougherty estimates for the effective yield stress and the effective consistency of the suspension: 



= r,V(l-<^)(l-^/(^„,)-2.5^™ (26) 

ry™ = r?(l-v.)^(l-^/^„)-^-^^("+^)^'" (27) 

Both the macroscopic yield stress and the macroscopic consistency diverge when Lp tends to- 
wards Lpra- The values obtained with Eq. (26) are plotted in Fig. 5, taking = 0.57. Good 
accordance is observed between theoretical and experimental data, accounting for the validity of 
Eq. (26). Note that no attempt to fit the value of has been performed. 

D Experimental results from the literature 

In this paragraph, the theoretical estimates are compared with experimental data of Ancey and 
Jorrot (2001), Erdogan (2005) and Geiker et al. (2002). 

The results for the yield stress are depicted in Fig. 7. The experimental data of Geiker et al. 
(2002) have not been drawn in Fig. 7 because they report to measure yield stress ten times larger 
than other researchers. 




Solid volume fraction 

Figure 7: Dimensionless yield stress r^°™/Tc vs beads volume fraction (/j for suspensions of Ancey 
and Jorrot (2001) (empty diamond-shaped) and Erdogan (2005) (solid diamond-shaped). The 
solid line is the theoretical prediction (26) with Lp^a = 0.635 and the straight dashed line is the 
dilute overall yield stress (24). The dashed curves are the theoretical Eq. (26) with = 0.67 
(dashed curve below the solid curve) and = 0.6 (dashed curve above the solid curve). 

The theoretical curve in Fig. 7 is calculated using = 0.635, as reported by Ancey and 
Jorrot (2001). We also plotted the curve using (pm = 0.67 (dashed curve below the solid curve) 
and ipm = 0.6 (dashed curve above the solid curve). The experimental data of Ancey and Jorrot 
(2001) are roughly well fitted by the theoretical estimate (26) with (pm = 0.635 even if the 
theoretical model always predicts yield stress lower than the measured ones. These discrepancies 
can not be impart to the experimental procedure. The suspensions was prepared by adding 
monodisperse glass beads to a clay dispersion and the yield stress was measured by means 
of a slump test, an experimental procedure which ensures that the materials remain isotropic 
and homogeneous. It is believed that the discrepancy between the experimental data of Ancey 
and Jorrot (2001) and our theoretical predictions could come from a wrong evaluation of the 
experimental maximum packing fraction. To support this opinion, it is worth noting that the 
estimate = 0.635 of Ancey and Jorrot (2001) was neither measured nor estimated from 
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experimental data. Moreover, this value 0.635 is greater than the values reported in the literature 
for a dense suspension of non colloidal suspension of monodisperse spherical particles [Ovarlez 
et al. (2006)]. It can be seen in Fig. 7 that the experimental data are well fitted by Eq. (26) with 

ifm = 0.6. 

Erdogan (2005) used a Couette-vane rheometer and applied a strong preshear to the samples 
before to measure the overall properties of the concrete suspension. Nevertheless, his results 
compare rather well to the theoretical ones, even if the lack of experimental data for values of 
the solid fraction larger than 0.4 does not allow to draw a definitve conclusion. 

We now pay attention to the results of the consistency measurements. Dimensionless con- 
sistency ■q^°'^/r] measured by Erdogan (2005) and Geiker et al. (2002) are plotted in Fig. 8 as 
a function of the normalized solid volume fraction (p/(pm- Four types of aggregates were used 
by Geiker et al. (2002) : glass beads, sea dredged, crushed, and a mix of 30% sea dredged and 
70% crushed aggregates. None of these coarse aggregates are monodisperse (the ratio of the 
bigger particle size to the smaller one is close to 4 for the all the aggregates). Moreover, only 
the glass beads are spherical. Then, most of the suspensions studied by Geiker et al. (2002) do 
not comply with the assumptions made to compute Eqs. (26) and (27). 
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Figure 8: Dimensionless consistency 7^1^°™ /r^ vs. the solid volume fraction ip for suspensions 
of Erdogan (2005) (solid diamond-shaped) and suspensions of Geiker et al. (2002) (empty circle 
for glass beads, other empty symbols for sea dredged, crushed and mixed aggregates). The solid 
line is the Krieger-Dougherty estimate for the consistency (n = 1) with = 0.65. The dashed 
lines are the same estimate with ipm = 0.55 (above the solid curve) and = 0.75 (below the 
solid curve). 

Consequently, it is not surprising that the data of Erdogan (2005) are relatively well fitted 
by the theoretical prediction while the fit between most of the experimental data of Geiker 
et al. (2002) and the theoretical estimate is less satisfactory. More precisely, while it seems that 
Eq. (27) is able to roughly estimate the consistency of suspensions of spherical particles (glass 
beads) even if they are not monodisperse, the discrepancy between experimental and theoretical 
data is no more acceptable for the other aggregates. As depicted in Fig. 8, tuning the value of 
the maximum solid fraction does not allow to improve the quality of the fit. 

It is believed that this discrepancy comes from the fact that the sea dredged, crushed and 
mixed aggregates are not spherical (an aspect ratio equal to 2 has been measured for all these 
particles). As the Krieger-Dougherty Eq. (23) applies only to suspensions of spherical particles, 
Eq. (27) cannot fit the properties of suspensions of non-spherical particles. One concludes that 
using a function g{(p) able to accurately estimate the linear properties of a suspension of non 
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spherical particles into Eqs. (18) and (19) will improve the quality of the estimates. This opinion 
is supported by the fact Geiker et al. (2002) observed that the relative consistency of non-spherical 
particles suspensions is greater than that of the glass beads suspensions when measured for the 
same solid volume fraction (see Fig. 8), a trend classically observed for Newtonian suspensions. 

V Validity of the theoretical model 

For all the results, the experimental data are more dispersed when the solid volume fraction 
increases and the discrepancy between the experimental and the theoretical results increases 
when the solid volume fraction tends towards that of the maximum packing fraction. It could 
be asserted that this discrepancy comes simply from the dispersion of the experimental results. 
Nevertheless, it is believed that heterogeneities of the secant modulus over the liquid domain are 
not any more negligible for the larger values of ^p. When the solid volume fraction value tends 
towards the maximum packing fraction, the fluid located near the closest points of two adjacent 
particles experiences a much larger shear rate than the fluid located far from these points (see 
Fig. 9). In this situation, it is necessary to take into account the heterogeneities of the secant 




Particles 
High shear rate 
Low shear rate 

Figure 9: Heterogeneities of the strain rate fleld for a concentrated suspension. 

modulus in the framework of more complex models in order to obtain more accurate estimates. 

To use an uniform estimate of the secant modulus induces a second drawback of the model. 
As no difference is made between all the points located in the fluid phase, such a model is unable 
to account for the fact that the fluid does not start to flow everywhere at the same time (or more 
precisely that the yield stress is not reached everywhere at the same time in the fluid domain). 

For the Bingham suspension, Castaheda (2003) has shown that beyond the yield stress, the 
behavior of the suspension is first nonlinear and then asymptotes to a purely linear behavior for 
higher strain rate values. As this phenomenon is induced by the presence of strong fluctuations 
of the strain rate fleld in the fluid phase at the onset of yield which induce fluctuations of the 
secant modulus, the simplifled model described above is not able to account for this nonlinearity. 
Nevertheless, it is possible to show that the simplifled approach performed in this paper furnishes 
useful estimates. 

Thanks to the properties of the minimization problem (8) it is easily established that the 
potential W (D) admits the following lower bound: 

W (D) > <^>i{D) + ^2{D) (28) 

with 



*.(D) = (l-^)..^min^(<i;^ (29) 
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and 

C&2 (£>) = (1 - (^) ^ min (30) 

The two quantities $1 and $2 define the eff^ective behavior of two suspensions with the same 
microstructure as the yield stress suspension. 

$1 defines the constitutive law of a suspension of particles embedded in a plastic material 
obeying the standard fiow rule for a von Mises criterion with yield shear stress Tc- ^2 is the 
overall potential of a suspension of rigid particles dispersed in a fiuid whose behavior is described 
by the incompressible viscous power-law: 

a = 2j](t-~^d-p6 trd = (31) 

Thanks to the properties of the function (^J™^ for ?n > 1 and of the set Q (D), it can be shown 
that the material defined by the potential is an isotropic rigid plastic material and the overall 
potential <^2 is that of an isotropic incompressible viscous material whose overall behavior is 
described by an homogeneous function of degree n. Moreover, if the suspending fluid obeys the 
Bingham law (n = 1), the effective potential $2 is that of a Newtonian fluid. 

Then, we can state that the macroscopic behavior of the isotropic suspension is energetically 
"underestimated" by a yield stress state law with same exponent as that of the bearing fluid. 
Both the overall yield criterion and the viscous behavior can be determined through the solution 
of two homogenization problems allowing to explicitly calculate the lower bound (28). Of course, 
the overall potential of the yield stress suspension is not the sum of two simple potentials. 
Nevertheless, it is obvious that the first term of the lower bound (28) prevails for the lower 
values of the macroscopic strain rate D, while for the larger values of the strain rate, the second 
term must be taken into account. This property allows to show that, although the macroscopic 
behavior of the suspension is not Herschel-Bulkley, it tends toward that of a yield stress fiuid 
with yield criterion defined by <l>i and viscous fiow rule described by ^2 for respectively the 
lower values of the macroscopic strain rate (near the yield stress) and the larger values of the 
macroscopic strain rate (far from the yield stress). These results are summarized in Fig. 10 




Figure 10: One-dimensional sketch of the functions W, and $2 as a function of the second 
invariant of D (left). The associated state laws (fictitious materials, lower bound and overall 
material) are depicted on the right part of the figure. The suspending fiuid is Bingham (n = 1). 



for a suspension of particles in a Bingham fluid. On the left part of Fig. 10 are drawn the 
overall potential W and both potentials $1 and $2 as a function of the second invariant of the 
macroscopic strain rate. For clarity, the lower bound ^1 + ^2 is not drawn. The overall behaviors 
of the suspension and of the materials associated with the potentials $1, ^2 and the lower bound 
$1 + $2 are depicted on the right part of Fig. 10. This diagram is very similar to the Fig. 3 of 
Castaheda (2003). S denotes the second invariant of the macroscopic deviatoric stress tensor. 
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Interestingly, using the method described in section III, it can be shown that Eqs (26) and 
(27) ahow to estimate the value of the yield stress defined by <I>i and the value of the consistency 
defined by <I>2. Then, even if the overall behavior of the suspension is not exactly Herschel- 
Bulkley, the estimates (26) and (27) can account for well defined properties of the suspension. 
Nevertheless it must be kept in mind that these estimates are not associated to lower bounds of 
the potentials $i and <I>2- A more accurate estimate for the overall behavior of a two-dimensional 
Bingham suspension, which is of the form suggested by the solid curve in the right part of Fig. 10, 
has been proposed by Castaheda (2003). 

It is clear from Fig. 10 that it is not possible to deduce the value of the macroscopic yield 
stress Sc from the measurements of the hight limit strain rate behavior of the overall material. 
Indeed, if one tries to evaluate the overall yield stress of the suspension from the high-limit strain 
rate behavior as depicted in Fig 10, the obtained estimate, denoted Ti^, is higher than the true 
yield stress Ec one must apply to make the suspension flowing. It seems that Geiker et al. (2002) 
used such a procedure to measure the overall yield stress of the concretes they studied. This 
result could explain why they measured yield stress values much higher (ten times) than the ones 
measured by others authors. 

The theory predicts that for the higher values of the macroscopic shear rate, the particles 
do not change the power-law index n. This result is in agreement with previous theoretical 
and experimental approaches to the rheological properties of concentrated suspensions in non- 
Newtonian fluid [Polinski et al. (1988); Castaheda (2003)]. 

VI Conclusion 

Theoretical estimates for the overall rheological properties of a suspension of noncolloidal and 
non-Brownian particles immersed in a nonlinear fluid have been proposed. To our knowledge, 
this approach and the associated estimates are original in the field of rheology. As the problem 
to be solved is nonlinear, it was necessary to make some approximations in order to compute 
simple analytical estimates. Here, the estimates are valid provided that the heterogeneities of 
the secant modulus can be neglected over the domain filled by the fluid phase. Then, the overall 
properties of the nonlinear suspension are estimated from that of a fictitious linear suspension 
having the same microstructure. 

In the case where the suspending fluid obeys a Herschel-Bulkley law, the estimates for the 
yield stress and the consistency read r^°™/rc = ^/ {1 — (p)g{ip) and i]^°"^/t] = g{ip)^'^^ / (1 — </?)"■" 
where g{f) denotes the ratio of the macroscopic to the microscopic properties of the fictitious 
linear material. 

Comparison of the theoretical results with experimental data for suspensions of monodisperse 
rigid noncolloidal particles dispersed in yield stress fluids have also been presented. The experi- 
mental procedure used to obtained these data was designed by Mahaut et al. (2007) to evaluate 
the purely mechanical contribution of the particles to the paste behavior, independently of the 
physicochemical properties of the materials. 

It was found from both the theoretical and the experimental approaches that the dimension- 
less yield stress and the dimensionless consistency depend on the bead volume fraction only. It 
has been shown that the yield stress/solid volume fraction and the consistency /solid volume frac- 
tion relationships are well fitted to the very simple laws r^°™/rc = (1 — ip)^^'^{l — ip/ipm)~^''^^'^"^ 
and ryhom/^ ^ _ ^)(i-«)/2 _ ^/^^)-i.25(n+i)c^„ respectively. 

As long as the coarse particles are spherical, the theoretical estimates agree quite well with 
experimental data found in the literature even if experimental procedures used by the authors do 
not exactly comply with the assumptions made to obtain the theoretical results. This is all the 
more satisfactory that only four numbers are necessary to calculate the estimates of the overall 
properties of the suspension as a function of the solid volume fraction: the close packing density 
of the particles and three suspending fluid properties (yield stress, consistency and power law 
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index) . 

We now plan to study the case of polydisperse systems and that of anisotropic particle 
distributions. To improve the theoretical estimates, it will also be necessary to separate the fluid 
domain in several geometrical domains in order to better describe the heterogeneities over the 
liquid domain. In particular this situation is expected to be founded for higher values of the 
solid volume fraction or for polydisperse suspensions. 
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